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SUMMARY 

In  this  study  the  effects  of  reduced  integration,  mesh  size,  and  element 
type  (i.e.  linear  or  quadratic)  on  the  accuracy  of  a  penalty  finite  element 
b,  ed  on  the  theory  governing  thick,  laminated,  anisotropic  composite  plates, 
are  investigated.  In  order  to  assess  the  accuracy  of  the  present  finite  ele¬ 
ment,  exact  closed-form  solutions  are  developed  for  cross-ply  and  antisymmetri 
angle-ply  rectangular  plates  simply  supported  and  subjected  to  sinusoidally 
distributed  mechanical  and/or  thermal  loadings,  and  free  vibration. 

INTRODUCTION 

The  increased  use  of  fi ber-reinforced  composites  in  aerospace  and  mech¬ 
anical  engineering  structures  is  primarily  due  to  their  high  stiffness-to- 
weight  ratio  and  due  to  their  anisotropic  material  property  that  can  be 
tailored  through  variation  of  the  fiber  orientation  and  stacking  sequence. 

With  the  increased  application,  the  thermo -mechanical  behavior  of  fiber- 
reinforced  composite  structures  is  receiving  greater  attention. 

Due  to  the  low  transverse  shear  moduli  relative  to  the  in-plane  Young's 
moduli,  the  transverse  shear  deformation  effects  are  more  pronounced  in 
plates  laminated  of  fiber-reinforced  materials  when  compared  to  thick  iso¬ 
tropic  plates.  A  shear  deformable  theory  that  is  an  extension  of  Reissner- 
Mindlin  theory  for  homogeneous  isotropic  plates  to  arbitrarily  laminated 
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anisotropic  plates  is  due  to  Yang,  Norris,  and  Stavsky  [1]  (see  also,  Whitney 
and  Pagano  [2]).  In  a  recent  study  the  author  [3,4]  derived,  using  the 
penalty-function  method,  the  Yang-Norris-Stavsky  (YNS)  theory  from  the  classical 
thin-plate  theory  (CPT)  of  laminated  plates.  A  finite-element  model  based 
on  the  YNS  theory  results  in  computationally  simpler  elements  compared  to 
previously  reported  elements  [5-3],  and  yet  possesses  competitive  accuracy. 

This  investigation  is  primarily  concerned  with  a  critical  study  of  the 
penalty  plate-bending  element  developed  by  the  author  [3,4].  Specifically, 
the  study  involves  an  investigation  of  the  numerical  accuracy  and  convergence 
of  the  element  when  applied  to  the  analysis  of  thick,  laminated  anisotropic 
composite  plates.  In  order  to  assess  the  accuracy  of  the  present  element, 
exact  closed-form  solutions  are  also  derived  for  freely  supported  rectangular 
plates  of  cross-ply  construction.  Comparison  is  also  made  with  other  finite- 
element  solutions,  and  3-D  elasticity  solutions  wherever  available.  Since 
the  penalty-function  method  is  used  in  developing  the  element,  the  so-called 
reduced  integration  [g  ,io]  comes  into  picture.  The  effect  of  reduced  inte¬ 
gration  on  the  accuracy  of  thin  and  thick  plates  is  also  investigated. 

SHEAR  DEFORMABLE  THEORY  OF  LAYERED  COMPOSITE  PLATES 


Consider  a  plate  of  uniform  thickness  t  composed  of  a  finite  number  of 
anisotropic  layers  with  arbitrary  orientations.  The  coordinate  system  is 
chosen  such  that  the  middle  plane,  F,  of  the  plate  coincides  with  the  x-y 
plane,  with  z-axis  normal  to  the  middle  plane. 

The  displacement  field  in  the  shear  deformable  theory  of  Whitney  and 

l  *1.1..  ^.c-ator 

Pagano  [2]  is  given  by 


u(x,y,z,T )  =  u0(x,y,x)  +  z^x(x,y,i) 
v(x,y,z,T)  =  vQ(x,y,T )  +  z^y(x,y,T) 
w  =  w(x ,y,t ) 
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where  u,  v,  and  w  are  the  displacements  along  x,  y  and  z  directions,  respec¬ 
tively,  uQ  and  vQ  are  the  in-plane  displacements  of  the  middle  plane,  and  ^ 
and  are  the  shear  rotations.  The  strain-displacement  equations  of  linear 
elasticity  become, 

ej  =  ex  =  UQ>X  +  »  e2  =  ey  =  v0,y  +  Z^y,y’ 

£c  =  Y  =  u  +v  +  z(t|i  +  U/  )  ,  (2 

6  xy  o,y  o,x  x,y  y,x; 

y  -ib  +  w,  ,  y  =  i i)  +  w,  ,  e  =  0 
rxz  vx  ’x  'yz  vy  y  z 

where  w,x  =  aw/ax t  etc. 

The  equilibrium  equations  of  the  theory  are  given  by. 
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where  R,  S,  and  I  are  the  normal,  coupled  normal -rot<.ry,  and  rotary  inertia 
coefficients, 

+/9  Zm+1 

f  W  ^  r  ^  m  \ 


(R,S,I)  =  [  (l,z,z2)  pdz  =  i  | 

J  -t/2  m  J  z 


(l,z,z2)  p^  dz 


p^m'  being  the  material  density  of  the  m-th  layer,  P:  and  P2  are  in-plane 
distributed  forces,  P3  is  the  transversely  distributed  force,  and  ,  and 

M.j  are  the  stress  and  moment  resultants  defined  by 

t/2  f t/2 

(N..M.)  =  (l,z)  o.  dz  ,  (Q.,Q2)  =  (axv,0vz^  dz 

1  1  -t/2  1  -t/2  xy  yz 
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Here  o..  (i  =  1,2,6)  denote  the  in-plane  stress  components  {ax  =  a x>  a2  =  ?y> 
and  as  =  axy) . 

Assuming  monoclinic  behavior  (i.e.,  existence  of  one  plane  of  elastic 
symmetry)  for  each  layer,  the  constitutive  equations  for  the  m-th  layer 
On  the  plate  coordinates)  are  given  by 
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so  that  (in  view  of  (5)  and  (6))  the  constitutive  equations  for  an  arbitrarily 
laminated  plate  are 
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The  plate  stiffnesses  A.-,  3..,  and  D.-  are  given  by 

l  J  1  J  I J 


m+1 


(Ai3*BiJ*Dij>  =  £ 


Q^O.z.z2)  dz  ,  (i,j  =  1,2,6) 
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W!"’  dz  ,  (=  -  6-1,  B  -  6-J;  i.j  =  4,5) 


where  are  the  stiffness  coefficients  of  the  m-th  layer  in  the  plate  coordi¬ 
nates,  and  zm  is  the  distance  from  the  mid-plane  to  the  lower  surface  of  the  m-th 
layer.  The  stress  and  moment  resultants,  wj  and  mT,  due  to  thermal  loading 
are  defined  by 


(n[,m{)  =  I 
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rZm+l 
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-  Qjj'  aJ(T0,zT1 )  dz  ,  (i,j  =  1,2,6) 


(9) 


where  a.  are  the  thermal  coefficients  of  expansion  in  the  plate  coordinates, 
and  T  is  the  temperature  chance  from  a  reference  state. 


T(x,y,z)  =  Tq(  x  ,y)  +  zT^x.y) 


(10) 


Note  that  the  temperature  variation  through  the  thickness  is  assumed  to  be 
linear,  consistent  with  the  plate  theory. 

Substituting  Eq.  (7)  into  Eq.  (3),  we  obtain  the  following  operator 
equation, 

[L ] { 6 }  =  {f }  (11) 


where  { o >  =  iu  ,v  ,w,»  [l]  is  the  (syrmetric)  matrix  of  differential 

oo  x  y 

operators. 


L11  =  Andn  +  2A16d12  +  A66d22  +  R  dTT 

Li 2  =  ^A12+A66^d12  +  A1 6dl 1  +  A26d22  ’  L13  =  0  ’ 

L14  =  Bndll  +  2S1 6dl 2  +  e66d22  +  S  hr 
L15  =  (B12+B66)d12  +  B1 6dl 1  +  B26d22  =  L24  ’ 

L22  =  2A26d12  +  A22d22  +  A66dll  +  R  hr  *  L23  =  0 

L25  =  2B26d12  +  S22d22  +  B66dll  +  S  dTx 
Ltj  =  ^dn1 -2A^cd-| 0-Accd0,  +  R  d  __  »  L,,,  =  -  A<cd, 
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(12) 


L35  =  _A45dl"A55d2’  L44  =  Dlldn+2D16d12+D66d22'A44  +  1  d 


TT 


L45  =  ^D12+D66^d12+D16dll+D26d22'A45’  L55  =  2D26dl 2+D22d22+D66dl 1  + 


and  the  components  of  the  generalized  force  vector,  { f } ,  are  given  by 

f ,  =  N  ‘  +  nI  +  Pi  ,  f,  =  n!  v  +  ill  +  p  „ 

1  1  ,x  6,y  1  2  6  ,x  z,y  r2 

f,  -P,  ,  f ,  =  m!  +  mJ  ,  f_  =  mI  +  M, 

J  3  4  l,x  6,y  5  6,x  2,y 

In  Eq.  (12),  d.  .  denote  the  differential  operators 

*  J 

dt*  3/s'-dfj  '  di  ■  “io  *  S’  •  ('.J  *  0,1,2) 


EXACT  CLOSED-FORM  SOLUTIONS 


The  boundary-val ue  problem  associated  with  the  equilibrium  of  layered 
anisotropic  composite  plates  involves  solving  the  operator  equation  (11) 
subjected  to  a  given  set  of  boundary  conditions.  It  is  not  possible  to 
construct  exact  solutions  to  Eq.  (11)  when  the  plate  is  of  arbitrary  geometry, 
constructed  of  arbitrarily-oriented  layers,  and  subjected  to  an  arbitrary  loading 
or  boundary  conditions.  However,  an  exact  closed-form  solution  to  Eq.  (11)  can 
be  constructed  when  the  plate  is  of  rectangular  geometry  with  the  following 
edge  conditions,  loading,  and  plate  construction  (cf.  [2]). 


Cross-Ply  Plates  in  Bendina  and  Vibration 


Boundary  conditions  (freely  supported) 
ujx.0)  =  ujx,b)  =  0,  N2(x,0)  =  N2(x,b) 

v</0.y)  =  vjla.y)  =  0,  N-|  (0,y)  =  N^a.y) 
w(x,0)  =  w(x,b)  =  w(0,y)  =  w(a,y)  =  0 


bx(x,0)  =  ii<x(x,b)  =  0,  M2(x,0)  =  M2(x,b)  =  0 
<l>y(0,y)  =  i|»y(a,yj  =  0,  M^O.y)  =  M^a.y)  =  0 


x 
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Load i nq  (sinusoidal) 

Pj3^  S1'n  ax  sin  gy,  T0=T^n  sin  ax  sin  gy,  T^=T^n  sin  a x  sin  gy 

_mn  (15) 

P1=P^n  sin  ax  cos  By  ,  P2=P^n  cos  ax  sin  By  ,  a  =  m*/a  ,  B  =  nn/b 

and  m  and  n  are  integers. 

Plate  construction  (cross-ply,  i.e.,  6m  should  be  either  0°  or  90°) 

A16  =  A26  =  A45  =  0  ’  B16  =  B26  =  0  ’  °16  =  °26  =  0  '  a6  =  0 

Linder  these  specific  conditions,  the  solution  (u  ,v  ,w,  ^ .  <P  )  to  Eq.  (11) 

oo  x  y 

is  of  the  form, 

uo  =  Umn  C0S  ax  sin  »  v0  =  vmn  sin  ax  cos  By 

w  =  ^mn  sin  ax  sin  By  (17) 


\  =  xmn  cos  ax  sin  Sy  ,  ^  =  Ymn  sin  ax  cos  By 


where  Umn,  Vmn>  etc.  are  parameters  to  be  determined  subjected  to  the  con¬ 


dition  that  the  solution  in  Eq.  (17)  satisfies  the  operator  equation  (11) 
Substituting  Eq.  (17)  into  Eq.  (11)  we  get 

[C](A}  =  (F } 


where 


{a}  =  (ij  v  w  x  v  (Ft  -  /pfT,n  f*710  pfnn  pmn.1 
mn5  mn*  mn’  Xmn’  Ymn}  •  {F}  "  {Fi  >F2  ’F3  »F4  >F5  5 


and  the  elements  of  the  coefficient  matrix,  [C],  are  given  by 
C11  =  'Alla  '  A66e2,  C12  =  ~^A1 2  +  A66^a6 


C1 3  ~  °’  C14 


4  =  -a2Bn  -  B66B2,  C15  =  -(B12  +  B66)ag 

2  a  2  r  ~  «  _  /> 


C22  =  “A22e  “  A66a  ’  C23  =  0>  C24  =  C15  ’ 

2  2  2  2 
C25  =  "B22e  ‘  A66a  *  C33  =  "a  A55  '  B  A44 

C34  =  "oA55’  C35  =  _eA44’  C44  =  'Dlla2  *  D66£2  *  A55 

C45  =  '(D12  +  °66)ae  *  C55  =  _D66a2  *  D22£2  '  A44 


.  A 
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■pfpn  _  Tflin  n.  .  V  =mu  ,  ,  n  \  =smi  T 

"l  “1  ^Allal  A1 2a2  To  ^Bnai  B12a2  T1  J 
=  Pf  +  S[(A12ttl  +  A22a2)  T^n  +  (Blzai  +  B^)  T?"  ] 

f4  =  at(Bnai  +  Bl2ct2^  'o  +  ^D11a1  +  °12a2^  = 

’  ^12*1  +  322a2)  C  +  (°12al  +  °22»2>  T  3 


=mn 


Thus,  for  a  given  a  =  m/s.,  s  =  mr/b,  generalized  forces  7<|,n  ,  and 
cross-ply  construction,  one  needs  to  solve  the  5  by  5  matrix  equation  (18) 
for  the  vector  {A}  of  amplitudes  of  the  generalized  displacements. 

Natural  vibration  frequencies  u  can  also  be  obtained  for  the  boundary 
conditions  and  cross-ply  construction  discussed  above  by  assuming  a  solution 
of  the  form  (cf.  Bert  and  Chen  [15]) 

uQ  =  U(x,y)  cos  ojt  ,  vQ  =  V(x,y)  cos  «  ,  w  =  W(x,y)  cos  bt 

(21) 

=  X(x,y)  cos  wt  ,  t|i  =  Y(x,y)  cos  bt 

a  y 


wherein  U(x,y),  V(x,y),  etc.  are  the  expressions  in  Eq.  (17).  Substituting 
Eq.  (21)  into  Eq.  (11)  (with  =  P2  =  Pi  =  0,  T  =  0) ,  we  get  the  eigenvalue 
equation, 

[C]  { A}  =  uj2[M]{a}  (22) 


where  the  elements  M.  .  =  M..  of  the  mass  matrix  [M]  are  given  by 

'  J  J  * 


M11  =  M22  =  M33  =  K  ’  M14  =  M25  =  3  ’  M44  =  M55  =  1 


(23) 


and  the  remaining  elements  are  zero. 


B.  Angle-Ply  Plates  in  Bending  and  Vibration 


Boundary  conditions  (simply  supported  edges) 
uQ( 0,y)  =  uQ(a,y)  =  0  ,  Ng(0,y)  =  Ng(a,y)  =  0 
vQ(x,0)  =  vQ(x,b)  =  0  ,  Ng(x,0)  =  Ng(x,b)  =  0 
w(x ,0)  =  w(x,b)  =  w(o,y)  =  w(a,y)  =  0 


(24) 


9 


( x . ° )  =  «x(x,b)  =  0  ,  M2(x,0)  =  M2(x,b)  =  0 

1>y(0,y)  *  *y(a,y)  =  0  ,  M^O.y)  ■  M^a.y)  =  0 


(24  cont. ) 


Loading  (sinusoidal) 

=  P^"  sin  ax  cos  By  ,  ~  ^2  cos  aX  sin  Sy  .  P3  =  sin  ax  sin  sy 

(25) 

x  -i-m  n  .  w  yfnn  *  . 

Tq  =  Tq  cos  ax  cos  ay  ,  T-j  =  T-j  sin  ax  sin  6y 
Plate  construction  (antisymmetric  angle-ply,  9/-e/s/-e,  ...) 

A16  =  A26  A45  =  0  '  Bn  =  B12  =  0  ’  °16  =  °26  =  0  ’  a6  =  0 


Then  the  solution  to  the  associated  boundary -value  problem  is  of  the  form 
(multiply  the  functions  with  cos  u;  for  eigenvalue  problem) 

u„  =  U  „  sin  ax  cos  Sy  ,  v„  =  V  cos  ax  sin  Sy 
0  mn  J  0  mn 

w  =  Wmn  sin  ax  sin  ay  (27) 

ti>  =  X  cos  ax  sin  By  ,  ti  =  Ym„  sin  ax  cos  By 
x  mn  y  mn 


where  a  =  mz/a  and  a  =  nz/b,  and  m  and  n  are  integers.  Substituting  Eq.  (27) 
into  Eq.  (11),  we  obtain, 


[K]  {a}  =  u2[M]  {a}  +  {F} 


(28) 


where 


'll 


{A} 

*66  1 


ill  „ , V  ,W  , V  , X  }  ,  and 
mn  ,m*  ,mn  mn  mn  ’ 


K,,  =  -(A^a2  +  A«*e2)  »  kt9  =  -<*6(A19  +  A, 


1 2 


66' 


^13  ~  B  =  "(*  +SBB?fi)/t,  K-|^  =  -2aBB-|fi/t 


'14 


>16 


26' 


1 5 


16' 


K22  =  A66  +  Z  A22^  ’  K23  =  0  ’  K24  =  -2aSB26/ t 

K25  =  -(*2b16  +  ^B25}/t  •  K33  =  *^55  +  s2a44} 

K34  =  “®A44 ’  ^35  =  _aA55^  ^  ’  ^44  =  B66  +  ®  B22  +  Aaa^^ 


44' 


'35 


>55' 


K45  =  -aS(D12  +  D66 )/-t  »  K55 


'44  ~  66  p  u22  1  44 

=  -(a20-||  +  S20fifi  +  A^/t2 


66 


55' 


The  elements  of  the  mass  matrix,  [M]  are  as  defined  in  Eq.  (23),  and  the 
elements  of  the  force  vector  are  given  by 


F,  =  -2  (Aucij  +  A12ot2)TQ  +  S ( B x 6ct j  +  B26a2)T; 
F2  =  -2  (A12oti  +  A22a2 ) Tq  +  ci ( B x. 6 ^ i  +  B2e:i2)Tj 


F3  =  P3/t 

Fs  =  -3  ( B i s ct x  +  B26a2)To  +  a(Duax  +  D12a2)Ti 
F4  =  -a  (Bieaj  +  B26a2)To  +  3(D32csj  +  D22a2)T-, 


(29) 


FINITE-ELEMENT  FORMULATION 


As  pointed  out  in  the  previous  section,  exact  solution  to  Eq.  (11) can 
be  obtained  only  under  special  conditions  of  geometry,  edge  conditions,  loadings, 
and  lamination.  Here  we  present  a  simple  finite-element  formulation  which  dees 
not  have  any  limitations  (except  for  those  implied  in  the  formulation  of  the 
governing  equations). 

Suppose  that  the  region  R  i  s subdi vided  into  a  finite  number  N  of  subregions 
or  finite  elements,  R  (e  =  1,2,...,N).  Over  each  element  the  generalized 
displacements  (u0>v0>w >tx>Vy)  are  interpolated  according  to 

:  i  :  i  s  2 

uo  -  :  Vi  •  v0  =  f  Vi  ’  w  =  t  Vi  ’ 


K  3  ¥  3 

z  *x1^  ,  -  Z  *y1*i 


(30) 


where  (a  =  1,2,3)  is  the  interpolation  function  corresponding  to  the  i-th  node, 
in  the  element.  Note  that  the  in-plane  displacements,  the  transverse  dis¬ 
placement,  and  the  slope  functions  are  approximated  by  different  sets  of 
interpolation  functions.  While  this  generality  is  included  in  the  formula¬ 
tion  (to  indicate  the  fact  that  such  independent  approximations  are  possible), 
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we  dispense  with  it  in  the  interest  of  simplicity  when  the  element  is  actually 

1  2  3 

programmed  and  take  =  ;■  (r  =  s  =  p).  Here  r,  s,  and  p  denote  the 

number  of  degrees  of  freedom  per  each  variable.  That  is,  the  total  number  of 
degrees  of  freedom  per  element  is  2r  +  s  +  2p. 

Substituting  Eq.  (30)  into  the  Galerkin  integrals  associated  with  tne 
operator  equation  (11^,  which  must  also  hold  in  each  element  Re> 


([L]{6>  -  {f}){*>  dx  dy  =  0 

J 13 


and  using  integration  by  parts  once  (to  distribute  the  differentiation 
equally  between  the  terms  in  each  expression),  we  obtain 


[K]  {v}  =  a)2  [M]  {7}  +  {FI¬ 


FO  r  static  bending,  Eq.  (29)  becomes 


[K11]  [K12]  [0]  [KU]  [K15]]  I  rui 


[K12]  [K22]  [0]  [K24]  [K25] 
[0]  [0]  [K33]  [K34]  [K35] 
[K14]  [K24]  [K34]  [K44]  [K45] 
[K15]  [K25]  [K35]  [K45]  [K55] 


where  the  {ul,  {v},  etc.  denote  the  columns  of  the  nodal  values  of  u,  v, 

~  ft 

respectively,  and  the  elements  (<*,S  =  1,2,. ..,5)  of  the  stiffness  matrix 
and  F.j  of  the  force  vector  are  given  by 

K!]  ■  A,lGlj  +  A,6<Glj  *  GJX()  *  A66G(j 
K!j2  ■  A12Glj  *  A16Gij  *  A26Glj  +  A66Gji 
‘I]  '  B,1HU  *  B16<HU  *  *  B66"(j 


V?j  +  Vlj  *  BZ6H<j  +  B66Hji 
VG)j  *  G**)  *  A22G^  ♦  A66G?j 

B16Hij  +  B66Ri j  *  B,2«^  *  »26H)j 
B26<Hlj  *  H*»)  *  B66H^  *  B22H^. 

A44SU  *  A45<^  *  Sj?)  *  AjjSfj 

A44R“  *  A45R1j  '  K?j  *  A45R“  *  ASSR)° 

DllT1j  +  °16<T1)  *  TjP  *  °66T1j  *  A44T?j 

D26(T1j  *  Tj(>  +  D66Tij  +  022I-'j  +  AssT^j 

:  R  j  9?  *■  dx  dy,  (a  =  £  =  1,2,3) 

Jr  1  J 
Ke 

1  s  {  ♦“  dx  dy  (a, 6  =  4,5) 

JR-  1  J 


?j  dx  dy  (ct  =  8  =4,5) 


f  $].  dx  dy,  (a 
a  ■  1 


=  1  ,2;  i  =  1 ,2,...,r) 


C 

\  2 

q?.  dx  dy,  (i  =  1,2,... ,s) 

J]R 

e 

fj,  dx  dy,  (a  =  4,5;  i  =  l,2,...,r) 

,  3  | 


where 


>!  *]. 
i,5  ^J,n 


dx  dy  ,  (i,j  =  1 ,2, . . . ,r) 


■j  3 

6.  a  .  dx  dy  ,  (i  s  l,2,...,r;  j  =  1,2,. 
>  *5  J  *n 


2  2 

.  dx  dy  ,  (1,j  =  1,2,...  ,s) 

i,5  J  ,n 
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Rf  •  -  !  9?  .  dx  dy  ,  (i  *  1,2 . s;  j  =  1,2 . t) 

ij  JR  i»5  J.ti 

Tj 1  =  [  dx  dy  ,  (i.j  =  1.2 . t)  ,  (?,n  =  O.x.y) 

ij  jp  i.s  J.n 

e 

xx  x  12  3 

and  G.  .  =  G . . ,  etc.  In  the  special  case  in  which  =  a.  =  all  of  the 

matrices  in  Eq.  (35)  will  coincide. 

In  the  present  study  rectangular  elements  with  four,  eight,  and  nine  nodes 
are  employed  with  the  same  interpolation  for  all  of  the  variables.  The 
resulting  stifness  matrices  are  20  by  20  for  the  A-ncde  element  and  40  by  40 
for  the  8-node  element.  As  pointed  out  in  a  recent  study  [6],  the  YNS 
theory  can  be  derived  from  the  corresponding  classical  thin-plate  theory 
by  treating  the  slope-displacement  relations 


3w  _  3w 

3x  "  ~"x  5  3y  "y 


as  constraints.  Indeed,  when  the  constraints  in  Eq.  (36)  are  incorporated 
into  the  classical  thin-plate  theory  by  means  of  the  penalty  method,  the 
resulting  equations  correspond  to  the  YNS  theory  with  the  correspondence, 


v  and  6  -  ± 

*x  y  *y 


It  is  now  well-known  that  whenever  the  penalty  method  is  used,  the  so-called 
reduced  integration  [  9]  must  be  used  to  evaluate  the  matrix  coefficients 


in  Eq.  (34).  That  is,  if  the  four-node  rectangular  element  is  used,  the  1  x  1 


Gauss  rule  must  be  used  in  place  of  the  standard  2x2  Gauss  rule  to  numerically 

evaluate  the  coefficients  K... 

■  J 


14 


NUMERICAL  RESULTS 


In  the  following,  we  present  numerical  results  for  bending  and  vibration 
analyses  of  laminated  anisotropic  plates.  In  most  of  the  problems  considered 
here  the  following  two  sets  of  material  properties  were  employed: 

Material  I:  Ej/Eg  =  25,  G12/E2  =  0.5,  G23/E2  =  0.2,  v]2  =  0.25 

(38) 

Material  II:  E-j/E 2  =  40,  G-j2/E2  =  0.6,  G23/E2  =  0.5,  =  0.25 

Further,  it  was  assumed  that  G-j2  =  G23  and  v^2  =  v-j3.  A  value  of  5/6  was  used 

2  2 

for  the  shear  correction  coefficients,  k.|  =  k^  (see  Whitney  [11]).  The  thermal 
properties  (i.e.,  coefficients  of  thermal  expansion,  cr )  are  given  during 
the  discussion  of  the  results.  All  of  the  computations  were  carried  on  an 
IBM  370/158  computer. 

Bending  Analysis 

Numerical  experiments  were  conducted  to  investigate  the  effect  of  element 
type  (i.e.,  linear,  eight -node  quadratic,  and  nine-node  quadratic),  mesh 
(L2  =  linear  2  by  2,  Q2  =  2Q8  =  8-node  quadratic  2  by  2  elements  in  quarter 
plate),  and  reduced  integration  on  the  deflections,  stresses,  and  natural  fre¬ 
quencies.  In  some  cases  the  difference  is  not  distinct  enough  to  show  on  the 
graph,  therefore  the  results  are  presented  in  tabular  form.  Table  1  shows  the 
effect  of  the  element  type,  mesh,  and  reduced  integration  on  the  accuracy  of 
the  maximum  deflection  and  stresses  (also  see  Fig.  1). 

In  Table  1  linear  and  quadratic  elements  are  compared  for  deflections 
and  stresses  of  a  three-layer  (t^=t3=t/4,  t2=t/2)  square  plate  (Material  I) 
subjected  to  sinusoidal  loading.  This  problem  is  also  equivalent  to  a  four- 
layer  (equal  thickness)  cross-ply.  The  stresses  were  computed  at  the  Gaussian 
points: 


M 


^inure  1  Closecl-forn  and  finite-el ement  solutions 
for  sir'pl y-suooorted  four-layer  cross-ply 
(Go/?no/903/f>°)  souare  plate  (Material  I) 
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w(f-,f,0)  E2tJ  x  10' 

u  =  _ — _ ~ _ 

5 


w  = 


V 


(7)  7-av(A,A,±t/2),  qo  =  P3 


'a'  q0  X' 


°"y  =  (a^  °y(A’A’±t/4)  *  V  (a)2  t  Txy(B’B’±t/2) 


(39) 


T’xz  =  txz(B’A)  in  layers  1  &  3;  T"yz  =  Tyz(A’B)  in  layer  2 


where  A  and  B  are  the  Gauss-point  coordinates  (as  a  fraction  of  sides  a 
and  b)  given  below: 

2x2  linear  4x4  linear  2x2  quadratic  4x4  quadratic 


A  0.125  0.0625  0.05283 

B  0.375  0.4375  0.4472 


0.02642 

0.4736 


In  Table  1,  R  denotes  reduced  integration,  F  =  full  integration,  and 
FR  =  full  integration  for  bending  terms,  and  reduced  integration  for  the 
shear  terms.  The  following  observations  can  be  made  from  the  results  of 
Table  1. 

1.  The  nine-node  element  gives  virtually  the  same  results  for  full 
(3x3  Gauss  rule),  reduced  (2x2  Gauss  rule),  and  mixed  (3x3  and 
2x2  Gauss  rules)  integrations.  However,  the  results  attained  by 
using  the  reduced  integration  are  the  closest  of  all  three  integra¬ 
tions  to  the  closed -form  solution  (CFS). 

2.  The  nine-node  element  and  the  eight-node  element  (with  reduced 
integration)  give  almost  identical  results  (contrary  to  the  belief 
that  the  nine-node  element  is  superior  to  the  eight-node  element) 
for  all  side-to-thickness  ratios.  In  fact,  the  eight-node  element 
is  relatively  more  accurate. 
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3.  Integration  rule  has  a  more  profound  effect  on  the  accuracy  in  the 
eight -node  element  than  in  the  nine-node  element.  Full  integration 
gives  less  accurate  results  than  the  reduced  integration,  end  the 
error  increases  with  the  side-to-thickness  ratio.  This  implies 
that  the  reduced  integration  is  a  must  for  thin  plates,  as  generally 
recognized  for  all  penalty  models. 

4.  Full  integration  results  in,  relatively,  smaller  errors  for  quad¬ 
ratic  elements  and  for  refined  meshes  than  for  linear  elements 
and/or  for  coarse  meshes.  That  is,  the  error  between  the  solutions 
obtained  by  2Q8-F  and  2Q8-R  is  smaller  than  those  obtained  by  4L-F 
and  4L-R,  and  the  error  between  the  solutions  obtained  by  4L-F  and 
and  4L-R  is  smaller  than  those  obtained  by  2L-F  and  2L-R. 

5.  Numerical  convergence  of  the  element  is  clear  from  the  results  (see 
also  Table  2  for  4Q8  results). 


In  the  results  to  be  presented,  reduced  integration  was  used  with  linear 
and  eight -node  quadratic  elements. 

In  Table  2,  3-D  elasticity  solution  (3-DES)  of  Pagano  and  Hatfield  [12],  the 
present  closed-form  solution  (CFS)  and  4Q8  finite -element  results  (FEM), 
and  the  finite -element  results  of  Panda  and  Natarajan  [8]  are  compared  for 
the  problem  discussed  in  Table  1.  The  present  finite-element  results  are  in 
excellent  agreement  with  the  closed-form  solution  of  the  equations  governing 
the  YNS  theory.  The  stresses  in  3-0  elasticity  solution  were  reported 
to  be  maximum  at  the  following  points: 


=  ax(|-,^,±t/2)  ,  0^  =  ay(^-,-£,±t/4)  , 


xy 

a 


Txy^.b^) 


Txz  =  Txz(a,I,0)  ’  Tyz  =  Tyz(2,b,0) 


(40) 
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Table  2  Comparison  of  maximum  deflection  and  stresses  in  three-ply 
(0°/90°/0°)  square  plate  subjected  to  sinusoidal  loading 
(t1=t3=t/4,  t0  =  t/2 ,  .Material  I,  4Q8-R  element) 


a/t 

1 

source 

{ 

w 

Gw 

- 

T  w 

- 

X 

y 

xy 

xz 

yz 

| 

3-DES  j 

51.189 

1.338 

0.835 

0.863 

0.153 

0.295 

2 

CFS 

50.625 

0.337 

0.725 

0.353 

0.1198* 

(0.2995) 

0.242 

1 

1 

FEM 

50.623 

0.337 

0.725 

0.353 

0.120 

(0.299) 

0.2227 

i 

i 

3-DES 

i 

19.537 

0.720 

0.663 

0.467 

0.21 9 

0.292 

4 

| 

CFS 

17.094 

0.403 

0.572 

0.306 

0.1388 

(0.347) 

0.195 

I 

I 

FEU 

17.094 

0.403 

0.572 

0.306 

0.1388 

(0.3472) 

0.195 

i 

l 

3-DES 

7.434 

0.559 

0.403 

0.276 

0.301 

0.196 

I  10 

CFS 

6.627 

0.496 

0.359 

0.240 

0.166 

(0.414) 

0.128 

j 

FEM 

6.627 

0.495 

0.359 

0.240 

0.165 

(0.414) 

0.128 

i 

Ref.  [8] 

6.299 

0.532 

0.307 

0.250 

— 

3-DES 

5.173 

0.543 

0.308 

0.230 

0.328 

0.156 

20 

CFS 

4.911 

0.524 

0.294 

0.219 

0.174 

(0.434) 

0.108 

FEM 

4 . 911 

0.524 

0.294 

0.219 

0.174 

(0.434) 

0.108 

Ref.  [8  ] 

4.847 

0.557 

0.307 

0.231 

— 

3-DES 

4.485 

0.539 

0.276 

0.216 

0.337 

0.141 

50 

CFS 

4.410 

0.533 

0.272 

0.213 

0.176 

(0.441) 

0.101 

FEM 

4.410 

0.533 

0.272 

0.213 

0.176 

(0.441) 

0.101 

Ref.  [8  ] 

4.421 

0.565 

0.287 

0.225 

— 

3-DES 

4.385 

0.539 

0.271 

0.214 

0.339 

0.139 

100 

CFS 

4.337 

0.535 

0.269 

0.212 

0.177 

(0.442) 

0.1001 

FEM 

4.336 

0.535 

0.269 

0.212 

0.177 

(0.442) 

0.1002 

Ref.  [8  ] 

4.363 

0.566 

0.284 

0.223 

-- 

— 

classical  thin-plate 
solution 

4.350 

0.539 

0.269 

0.213 

0.J39 

0.138 

_ 

*First  value  is  in  the  center  layer,  and  the  value  in  parenthesis  is  in  the 
top  and  bottom  layers. 


Note  that  the  present  finite -element  solutions  cannot  be  expected  to  agree 
well  with  the  3-D  elasticity  solution  since  the  element  is  based  on  the  YNS 
theory  rather  than  the  3-D  elasticity  equations.  The  accuracy  of  the  YNS 
theory  is  evident  from  the  results.  While  the  YNS  theory  seems  to  predict 
the  maximum  deflection  very  close  to  that  given  by  the  3-D  elasticity  solu¬ 
tion,  the  stresses,  especially  the  transverse  shear  stresses,  are  not  in 
good  agreement  for  thick  plates.  The  stresses  predicted  by  the  YNS  theory 
are  on  the  lower  side  of  the  3-D  elasticity  results. 

Table  3  shows  a  comparison  of  the  'exact'  solution  of  Pagano  [13], 
with  the  present  closed-form  solution  and  finite-element  results  for  a  three- 
layer  (equal  thickness)  cross -ply  square  plate  (Material  I).  Again,  the 
results  are  in  good  agreement  with  the  'exact'  solution.  Table  4  shows 
similar  results  for  a  rectangular  plate  (b/a  =  3)  of  the  same  construction 
as  above.  The  present  element  gives  more  accurate  results  compared  to  those 
obtained  by  Panda  and  Natarajan  [8  ],  who  used  a  complicated  shell  element 
in  their  analysis. 

Deflection  and  stresses  are  compared  in  Table  5  for  five-layer  (t-|=t3= 
t^=t/6,  t2=t^=t/4)  cross-ply  (0°/90o/0o/90°/0c’)  square  plate  (Material  I  ) 
under  sinusoidal  loading.  While  the  finite -element  solution  is  close  to 
the  closed-form  solution,  the  YNS  theory  gives,  relatively,  less  accurate 
solutions  for  the  five-layer  case  compared  to  the  three-layer  case  shown  in 
Table  2. 

Table  6  compares  the  maximum  deflection  and  stresses  in  a  three-layer 
(t-|=t3=0.1t,  t2=0.8t),  sandwich  square  plate  (0°/90°/0°)  subjected  to 
sinusoidal  loading.  The  material  of  the  face  sheets  (i.e.,  layers  1  and 
3) i s  the  same  as  Material  I,  and  the  core  material  is  transversely  isotropic 
with  respect  to  z  and  is  characterized  by  the  following  properties  (see 
Pagano  [13]), 


Table  4  Comparison  of  maximum  deflection  and  stresses  in  three-ply  (0°/90°/0c’) 
rectangular  (b=3a)  plate  subjected  to  sinusoidal  loading  (t.=t/3. 
Material  I,  4Q8-R) 


Source 


3-0  elasticity  [13] 


2  ! closed-form  sol 


present  FEM 


3-D  elasticity  [13] 


| closed -form 


present  FEM 


3-D  elasticity  [13] 


closed-form 


0.607 


0.0928 


0.0203  0.467 


0.0306 


Panda  &  Natarajan[8]|  0.752 


3-D  elasticity  [13]  |  0.610 


closed-form 


0.653 


0.65 


Panda  &  Natarajan[8]  0.565  0.654 


3-D  elasticity  [13]  |  0.520  0.628 


cl osed-form 


Panda  &  Natarajan[8]j  0.513  0.654 

_ _  I 


3-D  elasticity  [13]  |  0.508  0.62 


cl  osed-form 


0.506 


.506 


Panda  &  Natarajan  0.505  0.654 


classical  thin  plate  0>503  23 

solution 


36 


0.0299 


mm 


0.0287 


0.0259 


255 


0.0264 


253 


5 


0.0093 


0.00878 


0.0031 


0.0084 


0.00831 


0.00831 


0.0087 


0.00825 


0.00825 


0.0261  0.0086 


0.439 


0.0  19 


0.0134 


0.0134. 


1 


0.0110 


0.0108 


0.0  26 


0.0126 


0.0252  0.0083  0.440  0.0108 


Comparison  of  maximum  deflection  and  stresses  for  five-layer 
(0o/90o/0o/90o/0° )  square  plate  under  sinusoidal  loading 
(t,  =  t0  =  tc  =  t/6,  to  =  t,  =  t/4.  Material  I,  2Q8-R). 


a/t  i  source 


3-DES+[13]  53.415 


CFS 


.3902 


3-DES  [13]  18.668  0.685 


CFS  15.623  0.425 


5.602 


3-DES  [13] 


3  0. 


.53 


0.5355 


.633 


0.4879 


FEM  |  6.201  0. 


3-DES  [13] 


CFS 


3-DES  [13] 


.539 


3 


2 


.539 


.522  0.522 


.379 


3-DES  [13]  4.377  0.539 


.332  0 

.52' 

.3 

.52 

4.3505  0.539 


0.360 


.35 


0.359 


836 


235 


0.0236 


39 


22 


228 


0.0246 


0.0214 


222 


8 


0.0209 


0.0206 


0.0213 


6 


6 


0.0213 


.22 


.2 


0.275 


.238 


0.296 


0.258 


.2 


0.356 


.35 


0.272 


.35 


0.272 


0.102  i 


.229 


0.0932  i 


0.0692 
(0.173  )i 


0.205 


(0.172  | 


0.205 


+  3-DES  -  3-D  elasticity  solution 

*First  value  is  in  layers  1,3, and  5  and  the  value  in  parenthesis  is  in  layers 
2  and  4. 


100 


FEM 


3.065 


2.742 


^1^2  =  ^-0  »  ^3 /^2  =  12.5  (not  used  in  the  YNS  theory) 

(41) 

G1 3/G23  =  1-5E2  ’  G12  =  0,4  E2  ’  v12  =  v13  =  v23  =  0-25 

The  finite-element  results  are  in  good  agreement  with  the  closed-form  solu¬ 
tion;  however,  the  YNS  theory  seems  to  predict  the  stresses  quite  low  for 
thick  plates. 

A  comment  is  in  order  on  the  closed -form  solutions  for  the  cross -ply 
and  angle-ply  plates  subjected  to  sinusoidal  loadings  and  edge  conditions 
given  in  equations  (14)  and  (24).  In  the  case  of  cross-ply  plates  the 
assumed  solution  in  (17)  and  (21)  satisfies,  in  addition  to  the  boundary 
conditions  in  (14),  the  following  boundary  conditions , 

N-|  (x,0)  =  N-j  (x , b)  =  0;  N2(0,y)  =  N2(a,y)  =  0  (42) 

Similarly,  in  the  case  of  angle-ply  plates,  the  following  additional  boundary 
conditions  are  satisfied  by  the  solution  in  (27): 

M-j  (x,0)  =  M-j  (x , b)  =  0  ,  M2(0,y)  «  M2(a,y)  •  0  (43) 

From  the  variational  formulation  of  the  equations  in  (l)-(7),  it  follows  that 
the  only  boundary  conditions  that  are  physical  and  mathematically  correct 
are  those  in  (14)  and  (24).  That  is,  for  example,  if  u  is  specified  (essen¬ 
tial  boundary  condition)  at  a  boundary,  N-j  cannot  be  specified  (natural 
boundary  condition)  there,  and  vice  versa.  Similarly,  if  \p  is  specified 
at  a  boundary  then  M-|  cannot  be  specified  there.  Indeed,  in  the  finite- 
element  formulation  it  is  not  possible,  even  if  we  wish,  to  specify 
the  essential  and  natural  boundary  conditions  on  the  same  portion  of  the 
boundary.  As  a  result,  there  seems  to  exist  an  inherent  difference  between 
the  closed-form  solutions  and  finite -element  solutions.  However,  from  the 
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cross-ply  bending  results  it  is  clear  that  this  difference  is  not  noticeable. 

The  reason  could  be  attributed  to  the  fact  that  the  additional  boundary  condi¬ 
tions  (42)  are  also  nearly  met  by  the  finite  element  solution.  This  is  due 
to  the  fact  that  the  in-plane  displacements  are  zero  everywhere  in  the  plate, 
and  there  are  very  small  in-plane  forces  in  the  plate  (i.e.,  N-|  and  Mg  are  almost 
zero  everywhere  in  the  plate).  However,  for  thermal  loading,  N-j  and  Ng  are 
even  larger  (although  small),  and  the  results  may  not  agree  as  closely  as 
seen  thus  far. 

The  additional  boundary  conditions  (43),  in  the  case  of  angle-ply  plates, 
are  too  severe  to  be  satisfied  by  the  finite-element  method  in  all  cases. 

The  bending  moments  are  typically  large,  and  their  magnitude  decreases  with 
the  stiffness  coefficients,  B^g  and  B  .  These  coefficients  are  large  for 
an  even  number  of  layers  and  their  magnitude  decreases  with  increasing  number 
of  layers  (see  Tsai  [14]).  Thus,  in  angle-ply  plates  the  finite-element 
solutions  and  the  closed-form  solutions  could  differ  substant’’al ly  for  small 
number  of  layers  (say  n  =  2)  and  be  almost  equal  for  large  number  of  layers 
(n  -  y>) . 

One  other  comment  that  applies  to  the  closed-form  solution  (27)  of  angle- 
ply  plates  is  that  the  assumed  solution  forces  u  and  v  to  vanish  at  y  =  b/2 
and  x  =  a/2,  respecti vely: 

u(x,b/2)  =  0  ,  v  (a/2,y)  =  0  (44) 

In  the  finite -element  modeling  of  a  quarter  plate,  the  following  symmetry 
conditions  were  used  for  all  edge  conditions  and  laminations: 

u(a/2,y)  =  0,  y^(a/2,y)  *  0,  v(x,b/2)  =  0,  iy(x,b/2)  =  0  (45) 
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Clearly,  the  symmetry  boundary  conditions  in  (45)  are  physical  whereas 
those  in  (44)  are  not  physically  meaningful.  This  marks  another  difference 
(which  turns  out  to  be  a  major  one)  between  the  closed-form  and  finite-element 
formulations.  Of  course,  one  can  force  the  symmetry  boundary  conditions  (44) 
in  the  finite-element  model.  Indeed,  when  the  symmetry  conditions  in  (44) 
are  forced;  the  finite-element  solution  satisfies  the  boundary  conditions  in 
(43),  and  consequently  the  two  solutions  agree  very  closely. 

The  above  observations  are  supported  by  the  numerical  results  obtained 
for  the  angle-ply  plates.  Table  7  shows  the  maximum  nondimensionalized  de¬ 
flection  for  two-  and  sixteen-layer  angle-ply  (e/-e/e/-e , . . . )  square  plate 
(Material  II)  under  sinusoidal  loading.  Since  the  stiffness  coefficient 
B-J6  (-1/n)  is  the  largest  at  e  =  45°  for  a  given  even  number  of  layers,  it 
is  the  worst  case  in  terms  of  the  agreement  with  the  closed -form  solution. 
However,  for  n  =  16,  the  finite-element  solution  is  in  good  agreement  with 
the  closed-form  solution.  As  expected,  the  results  for  e  =  5°  and  30°  are 
closer  to  the  closed-form  solution,  even  for  n  =  2,  compared  to  the  results 
obtained  for  9  =  45°. 

Figure  2  shows  the  effect  of  side-to-thickness  ratio  (a/t),  number  of 
layers,  lamination  angle,  and  the  symmetry  boundary  conditions  on  the  non¬ 
dimensionalized  maximum  deflection.  Note  that  for  e  =  45°,  the  finite -element 
solution  obtained  with  the  physical  boundary  conditions  in  (45)  is  about  one- 
half  of  the  closed-form  solution  for  n  =  2,  whereas  for  n  =  16  they  are 
almost  identical.  When  the  symmetry  boundary  conditions  in  (44)  are  employed, 
instead  of  those  in  (45),  in  the  finite-element  model,  the  closed -form  and 
finite-element  solutions  are  in  excellent  agreement  for  any  n  and  e  (Figure 
2  contains  results  only  for  n  =  4  and  e  =  45°). 


Si 


a/t  10  20  30  40  50 


Figure  2  Effect  of  side-to-thickness  ratio,  number  of  layers, 
lamination  angle,  and  the  symmetry  boundary  conditions 
on  the  nondimensional  ized  deflection  of  an  angle-ply 
sauare  nlate  (Material  II). 
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Figure  3  shows  the  nondimens ional ized  deflection  versus  the  lamination 
angle  for  two,  four  and  sixteen  layer  angle-ply  (0/-0/0/-0/. . ./)  square  plates 
(Material  II)  under  sinusoidal  loading.  Dark  lines  and  symbols  correspond, 
respectively,  to  the  closed-form  solution  and  the  finite -element  solution 
with  the  symmetry  conditions  of  the  CFS  imposed.  They  are  in  excellent  agree¬ 
ment  with  each  other.  The  broken  lines  in  Figure  3  correspond  to  the  finite- 
element  solution  obtained  by  using  the  physically  appropriate  symmetry  boundary 
conditions.  Open  symbols  are  used  for  thick  plates  (a/t  =  10)  and  dark  symbols 
are  used  for  thin  (a/t  =  100)  plates. 

Table  8  shows  the  nondimensionalized  deflection  as  a  function  of  side- 
to-thickness  ratio,  element  type  and  mesh  for  a  cross-ply  (0°/90°/0°,  t.=t/3) 
plate  (Material  I  with  ^  =  3a-|)  under  sinusoidal  temperature  distribution. 

The  finite-element  results  are  in  excellent  agreement  with  the  closed -form 
solution.  Figure  4  shows  the  effect  of  side-to-thickness  ratio  (a/t) 
on  the  nondimensionalized  deflection  and  stresses  for  cross-ply  and  angle-ply 
plates  subjected  to  thermal  loading  (P^  =  0).  The  scale  for  the  deflection 
and  stresses  is  amplified  to  show  the  effect  of  thickness-shear  strain. 

Free  Vibration  Analysis 

Once  again,  the  effect  of  reduced  integration  and  the  use  of  eight  and 
nine  node  elements  on  the  accuracy  of  the  natural  frequencies  are  studied 
using  three-layer  cross-ply  (0°/90°/0°)  example  problem  of  Table  1.  Table  9 
shows  the  nondimensionalized  fundamental  natural  frequency  as  a  function  of 
side-to-thickness  ratio,  type  of  integration,  and  type  of  element.  From  the 
results  obtained,  it  is  clear  that  no  single  integration  type  is  the  best 
choice  (in  contrast  to  the  observations  made  in  the  bending  analysis)  for 
all  ratios  of  side  to  thickness.  However,  the  reduced  integration  still 
seems  to  be  giving  better  results  for  most  ratios  of  side  to  thickness. 


e  10°  203  30°  40°  50° 

Figure  3.  Nondimensional i zed  deflection  versus  the 


lamination  angle  (e)  for  antisymmetric  angle- 
ply  square  plates  (Material  II) 
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Table  8  Three-Layer  (0°/90°/0°)  Simply-Supported  Square  Plate 

(Material  I ,  ct2  =  3ct-| )  Subjected  to  Sinusoidal  Temperature 


Distri bution 


a/b=l/3 


L2x2 

1.1154 

L4x4 

1 .0730 

Q2x2 

1.0617 

Q4x4 

1.0592 

1 

1.0593 

1.1158 

L4x4 

1 .0733 

Q2x2 

1 .0604 

Deflection 


a/b=0.5 


1.1311 
1.0380 
1.0727 
1 .0740 
T.0741 


a/b=l  .0  a/b=l  .5  a/b=2.0 


1.1530 
1 .1090 
1.0917 
1.0948 
T.0949 


131 

9  1 

088 

7  1 

074 

0  1 

074 

7  1 

1 .0596 


1.0748 


1.0963 


1.0369 

0.99741 

0.98154 

0.98458 

0.9847 


1.0374 

0.99790 

0.98358 

0.98506 


0.9851 


0.80483 

0.77422 

0.76206 

0.76427 

0.7643 


0.80416  ' 
0.77362  i 
0.76259  | 
0.76370 


0.7638 


L2x2 

1 .1 

18 

6  1  .1 

374 

1 .1 

658 

L4x4 

75 

8  1.0 

936 

1.1 

204 

02x2 

1.0 

62 

789 

1.1047 

1.0411 

0.79992 

1.0012 

0.76982 

0.98717 

0.75924 

0.98821 

0.76003 

1.0619 


1.1281 
1 .0842 
1 .0704 
1.0700 
1.0701 


1.1625 
1 .1150 
1 .1001 
1.0998 


1.0795 


1.1565 
1.1106 
1.0953 
1 .0958 
1 .0959 


1.2228 
1.1703 
1 .1529 
1.1534 


1 .1058 


1 .2012 
1 .1522 
1 .1354 

1 .1364 

1.1365 


1.2972 
1.2406 
1.2213 
1 .2224 


0.9883 


1.0513 
1 .0104 
0.99620 
0.99718 
0.9973 


1 .0713 
1.0293 
1.0147 
1 .0157 


0.7601 


0.78946  j 
0.76022  ! 
0.74996  j 
0.75070  j 
0.7508  | 


0.77331  j 
0.74474 
0.73473 
0.73545  ' 


( 0 ° / 90 0 / 90 0 / 0 0 ) ,  Material  II 


(457-457457-45°),  Material  II 

a 


(45  7-45  745  7-45°) ,  Material  I 

b  =  2 
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(b)  stresses  (at  Gauss  points) 
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Figure 4  Cffect  of  thickness  on  the  thermal  response  of 
cross-ply  and  angle-ply  simply  supported  plates 
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Figure  5a  shows  the  effect  of  side-to-thickness  ratio  on  the  nondimen- 
sionalized  fundamental  frequency  of  cross-ply  rectangular  plates  (Material  II). 
Similar  results  are  presented  for  angle-ply  plates  in  Figure  5b.  Note  again, 
that  the  frequencies  predicted  by  the  physical  symmetry  boundary  conditions 
(in  FEM)  are  substantially  different  from  the  closed-form  solution  for  two- 
layer  (45°/-45°)  angle-ply  plates.  However,  for  eight-layer  angle-ply  plates, 
this  difference  vanishes.  Figure  5b  also  shows  the  nondimensionalized  frequency 
versus  the  lamination  angle  for  four-layer  angle-ply  plates  (a/b  =  1,  a/t  =  10, 
Material  II).  The  finite-element  solutions  with  both  types  of  symmetry 
boundary  conditions  are  close  to  the  closed-form  solution  and  the  differences 
can  be  seen  on  the  plot. 


SUMMARY  AND  CONCLUSIONS 

In  this  study  we  investigated  the  effect  of  mesh,  element  type,  numerical 
integration  (full  and  reduced),  and  boundary  conditions  on  the  accuracy  of 
deflections,  stresses  and  natural  frequencies  associated  with  cross-ply  as 
well  as  antisymmetric  angle-ply  plates  by  comparing  the  solutions  with  the 
closed-form  solutions  developed  herein.  The  finite-element  solutions  are 
found  to  be  in  close  agreement  with  the  closed-form  solutions  for  only  2  by  2 
mesh  of  quadratic  elements  in  the  quarter  plate  (for  4  by  4  mesh,  they  are 
almost  identical  to  the  corresponding  closed-form  solutions).  The  closed- 
form  solution  for  angle-ply  plates  imply  nonphysical  symmetry  boundary  condi¬ 
tions,  and  its  effect  on  deflections  and  natural  frequencies  is  discussed. 
Reduced  integration  is  essential  for  the  analysis  of  thin  plates,  but  it  is 
not  crucial  for  thick  plates. 

The  element  developed  herein  has  been  applied  to  the  analysis  of  bi- 
modulus  (materials  having  different  properties  in  tension  and  compression) 


Fioure  5  Effect’ of  side-to-t'rickness  ratio  on  the  non- 

dimensional  ized  natural  frequencies  of  composite 
plates  (Material  II) 
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composite  plates  by  Bert  et  al.  [16],  and  Reddy  and  Chao  [17].  Application 
of  the  present  element  to  transient  dynamic  analysis  is  awaiting. 
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